################################################################################
### Politicians' Complaint Response: E-Governance and Personal Relationships ###
### William O'Brochta                                                        ###
### Louisiana Tech University                                                ###
################################################################################
#This file replicates all figures and tables in the main text and the
#supplemental information except: Figures 4.1-4.3 and Tables 4.4-4.5.
#See MC Question Text Analysis.R to replicate those figures and tables.


normalize<-function(x){(x-min(x, na.rm=T))/(max(x, na.rm=T)-min(x, na.rm=T))}

library(ggplot2)
library(plm)
library(fastDummies)
library(plyr)
library(dplyr)
library(tidyr)
library(car)
library(lmtest)
library(sandwich)
library(multiwayvcov)
library(lmerTest)
library(rstanarm)

#Read and merge Delhi
praja<-read.csv("MCD_Praja.csv", header=T, stringsAsFactors = F)
constituency<-read.csv("MCD_ConstituencyData.csv", header=T, stringsAsFactors = F)
survey<-read.csv("MCD_2019_Survey.csv", header=T, stringsAsFactors = F)

prajaconstituency<-merge(praja, constituency, by="Constituency", all.x=T)
combined<-merge(prajaconstituency, survey, by=c("Constituency", "Year"), all.x=T)

combined$female<-ifelse(combined$Gender=="Female",1,0)
combined$BJP<-ifelse(combined$Party=="BJP",1,0)
combined$AAP<-ifelse(combined$Party=="AAP",1,0)
combined$SS<-ifelse(combined$Party=="SS",1,0)
combined$SC_Pct<-combined$SC_Population/combined$Population
combined$MOV<-combined$MOV_Votes/(combined$WinningVotes+combined$RunnerUp_Votes)
combined$BA<-ifelse(combined$Education=="BA",1,0)
combined$Above10<-ifelse(combined$Education!='Below10',1,0)

#Predicting score based on demographics
delhi2019c<-combined[combined$City=='Delhi' & combined$Year=='2019',]
delhi2018c<-combined[combined$City=='Delhi' & combined$Year=='2018',]

delhi2019c$TotalScore_n<-normalize(delhi2019c$TotalScore)
delhi2019c$Attendance_n<-normalize(delhi2019c$Attendance)
delhi2019c$IssuesRaised_n<-normalize(delhi2019c$IssuesRaised)
delhi2019c$AwarenessAccessibility_n<-normalize(delhi2019c$AwarenessAccessibility)
delhi2019c$PerceivedPerformance_n<-normalize(delhi2019c$PerceivedPerformance)

delhi2018c$TotalScore_n<-normalize(delhi2018c$TotalScore)
delhi2018c$Attendance_n<-normalize(delhi2018c$Attendance)
delhi2018c$IssuesRaised_n<-normalize(delhi2018c$IssuesRaised)
delhi2018c$AwarenessAccessibility_n<-normalize(delhi2018c$AwarenessAccessibility)
delhi2018c$PerceivedPerformance_n<-normalize(delhi2018c$PerceivedPerformance)


modelc1<-lm(TotalScore_n~female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+MOV+
             Ward
             , data=delhi2019c)
coefc1<-coeftest(modelc1, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2019c$Ward, type="wild"))

modelc1.1<-lm(TotalScore_n~female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+MOV+
             Ward
           , data=delhi2018c)
coefc1.1<-coeftest(modelc1.1, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2018c$Ward, type="wild"))


modelc1.2<-lm(Attendance_n~female+SC_Pct+BJP+
              ReservedCaste+Population+C_Standing+
              NumberCommittees+Age+BA+MOV+
              Ward
            , data=delhi2019c)
coefc1.2<-coeftest(modelc1.2, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2019c$Ward, type="wild"))

modelc1.3<-lm(Attendance_n~female+SC_Pct+BJP+
                ReservedCaste+Population+C_Standing+
                NumberCommittees+Age+BA+MOV+
                Ward
              , data=delhi2018c)
coefc1.3<-coeftest(modelc1.3, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2018c$Ward, type="wild"))



modelc1.4<-lm(IssuesRaised_n~female+SC_Pct+BJP+
              ReservedCaste+Population+C_Standing+
              NumberCommittees+Age+BA+MOV+
              Ward
            , data=delhi2019c)
coefc1.4<-coeftest(modelc1.4, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2019c$Ward, type="wild"))

modelc1.5<-lm(IssuesRaised_n~female+SC_Pct+BJP+
                ReservedCaste+Population+C_Standing+
                NumberCommittees+Age+BA+MOV+
                Ward
              , data=delhi2018c)
coefc1.5<-coeftest(modelc1.5, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2018c$Ward, type="wild"))



modelc1.6<-lm(AwarenessAccessibility_n~female+SC_Pct+BJP+
              ReservedCaste+Population+C_Standing+
              NumberCommittees+Age+BA+MOV+
              Ward
            , data=delhi2019c)
coefc1.6<-coeftest(modelc1.6, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2019c$Ward, type="wild"))

modelc1.7<-lm(AwarenessAccessibility_n~female+SC_Pct+BJP+
                ReservedCaste+Population+C_Standing+
                NumberCommittees+Age+BA+MOV+
                Ward
              , data=delhi2018c)
coefc1.7<-coeftest(modelc1.7, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2018c$Ward, type="wild"))



modelc1.8<-lm(PerceivedPerformance_n~female+SC_Pct+BJP+
              ReservedCaste+Population+C_Standing+
              NumberCommittees+Age+BA+MOV+
              Ward
            , data=delhi2019c)
coefc1.8<-coeftest(modelc1.8, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2019c$Ward, type="wild"))

modelc1.9<-lm(PerceivedPerformance_n~female+SC_Pct+BJP+
                ReservedCaste+Population+C_Standing+
                NumberCommittees+Age+BA+MOV+
                Ward
              , data=delhi2018c)
coefc1.9<-coeftest(modelc1.9, vcov=function(x) vcovBS(x, R=5000, cluster=delhi2018c$Ward, type="wild"))

#Table SI.2.1
library(stargazer)
stargazer(modelc1, modelc1.1, modelc1.2, modelc1.3, modelc1.4,
          modelc1.5, modelc1.6, modelc1.7, modelc1.8, modelc1.9,
          se=list(coefc1[,2], coefc1.1[,2], coefc1.2[,2],
                  coefc1.3[,2], coefc1.4[,2], coefc1.5[,2],
                  coefc1.6[,2], coefc1.7[,2], coefc1.8[,2], coefc1.9[,2]), digits=3)


#Questions and Complaints
delhiquestions<-read.csv("MCD_Questions.csv", header=T, stringsAsFactors = F)
delhicomplaints<-read.csv("MCD_Complaints.csv", header=T, stringsAsFactors = F)
delhiquestions$Date<-as.Date(delhiquestions$Asked_Date)
delhicomplaints$Date<-as.Date(delhicomplaints$Creation_Date)


#Aggregate questions and complaints at ward level
#Ward committe membership doesn't change by year; 
#questions/complaints run 2018-04-01 to 2019-04-01 so use 2019 councillor report card

delhiquestions<-delhiquestions[delhiquestions$Date>"2018-04-01",]
delhiquestions<-delhiquestions[delhiquestions$Date<"2019-04-01",]
delhicomplaints<-delhicomplaints[delhicomplaints$Date>"2018-04-01",]
delhicomplaints<-delhicomplaints[delhicomplaints$Date<"2019-04-01",]

#Take only questions asked in ward committees that are ward related
delhiquestions<-dummy_columns(delhiquestions, select_columns="Committee_Name")
delhiquestions$Ward_Committee<-rowSums(delhiquestions[,c(17,18,24,25,28:31,33:35,39)])

delhiwardquestions<-delhiquestions[delhiquestions$Ward_Committee==1,]
delhiwardquestions<-delhiwardquestions[delhiwardquestions$Ward_Related=='Yes',]

complaintscoding<-read.csv("ComplaintsCoding.csv", header = T, stringsAsFactors = F)
questionscoding<-read.csv("QuestionsCoding.csv", header=T, stringsAsFactors = F)
reducedcoding<-read.csv("ComplaintsReducedCoding.csv", header=T, stringsAsFactors = F)

#Merge by coding for questions and complaints
delhiquestions2<-merge(delhiquestions, questionscoding, by="Level3", all.x=T)
delhiquestions3<-merge(delhiquestions2, reducedcoding, by="ComplaintCoding", all.x=T)

delhiwardquestions2<-merge(delhiwardquestions, questionscoding, by="Level3", all.x=T)
delhiwardquestions3<-merge(delhiwardquestions2, reducedcoding, by="ComplaintCoding", all.x=T)
delhicomplaints2<-merge(delhicomplaints, complaintscoding, by="Complaint_Type", all.x=T)
delhicomplaints3<-merge(delhicomplaints2, reducedcoding, by="ComplaintCoding", all.x=T)

delhiquestions3<-dummy_columns(delhiquestions3, select_columns="Category")
delhiwardquestions3<-dummy_columns(delhiwardquestions3, select_columns="Category")
delhicomplaints3<-dummy_columns(delhicomplaints3, select_columns="Category")


#Analyze the Questions
#save(delhiquestions3, file="delhiquestions3.RData")

#Now aggregate across category for each ward
delhiwardquestions4<-delhiwardquestions3[,c(5,45:61)]
delhicomplaints4<-delhicomplaints3[,c(3, 10:24)]

delhiwardquestions5<-delhiwardquestions4%>%
  group_by(Committee_Name)%>%
  summarise_all(sum, na.rm=T)
delhiwardquestions5<-delhiwardquestions5[-13,]

#Table SI.4.3
library(xtable)
print(xtable(delhiwardquestions5, type='latex', row.names=F), file='delhiwardquestions5.tex', include.rownames=FALSE)


#Table SI.3.5
delhicomplaints5<-delhicomplaints4%>%
  group_by(Ward)%>%
  summarise_all(sum, na.rm=T)
delhicomplaints5<-delhicomplaints5[-13,]

delhicomplaints5_pct<-delhicomplaints5
delhicomplaints5_pct[,2:16]<-delhicomplaints5_pct[,2:16]/rowSums(delhicomplaints5_pct[,2:16])

print(xtable(delhicomplaints5_pct, type='latex', row.names=F), digits="2", include.rownames=FALSE)



delhiwardquestions6<-gather(delhiwardquestions5, category, count, 
                            Category_Drainage:Category_NA, factor_key=TRUE)

delhicomplaints6<-gather(delhicomplaints5, category, count, 
                         Category_Drainage:Category_NA, factor_key=TRUE)

delhicombined<-merge(delhiwardquestions6, delhicomplaints6, 
                     by.x=c("Committee_Name", "category"),
                     by.y=c("Ward", "category"), all.x=T, all.y=T)

delhicombined[is.na(delhicombined)] <- 0 
colnames(delhicombined)<-c("Ward", "Category", "Questions", "Complaints")

delhicombined2<-delhicombined[,-2]%>%
  group_by(Ward)%>%
  summarise_all(sum, na.rm=T)

delhicategory<-delhicombined[,-1]%>%
  group_by(Category)%>%
  summarise_all(sum, na.rm=T)









#Aggregate by corporator
personquestions<-as.data.frame(table(delhiquestions$Constituency))
personquestions2<-as.data.frame(table(delhiwardquestions$Constituency))
personquestions3<-merge(personquestions, personquestions2, by="Var1", all.x=T, all.y=T)
colnames(personquestions3)<-c("Constituency", "Questions", "WardQuestions")
personquestions3<-personquestions3[order(-personquestions3$Questions),]
personquestions3$WardQuestions<-ifelse(is.na(personquestions3$WardQuestions),0,personquestions3$WardQuestions)

#Table SI.4.1
personquestions4<-personquestions3[1:25,]
print(xtable(personquestions4, type='latex', row.names=F, digits=0), include.rownames=FALSE)

#Table SI.4.2
personquestions3<-personquestions3[order(-personquestions3$WardQuestions),]
personquestions4<-personquestions3[1:25,]
print(xtable(personquestions4, type='latex', row.names=F, digits=0), include.rownames=FALSE)


#Ward Committee-Issue Level
delhicombined$Questions_Pct<-delhicombined$Questions/sum(delhicombined$Questions)
delhicombined$Complaints_Pct<-delhicombined$Complaints/sum(delhicombined$Complaints)
delhicombined$Questions_n<-normalize(delhicombined$Questions)
delhicombined$Complaints_n<-normalize(delhicombined$Complaints)
delhicombined$logComplaints<-log(delhicombined$Complaints+0.001)
delhicombined$logQuestions<-log(delhicombined$Questions+0.001)
delhicombined$Complaints_Question<-delhicombined$Complaints/delhicombined$Questions

delhicategory$Complaints_Question<-delhicategory$Complaints/delhicategory$Questions





###Check Ward Classification
wardlist<-read.csv("Constituency_Ward.csv", header=T, stringsAsFactors = F)

osm<-read.csv("delhicomplaints_rafael_openstreetmap_01.csv", header=T, stringsAsFactors = F)
osm$Constituency<-ifelse(osm$Constituency=="not found", NA, osm$Constituency)
osm$Constituency2<-paste0(substr(osm$Constituency, start = 1, stop = 3),"-")
osm$Constituency2<-paste0(osm$Constituency2,substr(osm$Constituency, start = 4, stop = 4))
osm$Constituency2<-ifelse(osm$Constituency2=="NA-NA",NA,osm$Constituency2)
osm<-osm[,-4]
osm2<-merge(osm, wardlist, by.x="Constituency2", by.y="Constituency", all.x=T, all.y=F)
colnames(osm2)<-c("OSMConstituency", "Ward", "Complaint_No", "Area", "OSMConfidence", "OSMWard")


google<-read.csv("delhicomplaints_rafael_google.csv", header=T, stringsAsFactors = F)
google$Constituency<-ifelse(google$Constituency=="not found" | 
                            google$Constituency=="outside delhi", NA, google$Constituency)
google$Constituency2<-paste0(substr(google$Constituency, start = 1, stop = 3),"-")
google$Constituency2<-paste0(google$Constituency2,substr(google$Constituency, start = 4, stop = 4))
google$Constituency2<-ifelse(google$Constituency2=="NA-NA",NA,google$Constituency2)
google<-google[,-4]
google2<-merge(google, wardlist, by.x="Constituency2", by.y="Constituency", all.x=T, all.y=F)
colnames(google2)<-c("GoogleConstituency", "Ward", "Complaint_No", "Area", "GoogleConfidence", "GoogleWard")


osmgoogle<-merge(osm2, google2, by=c("Ward", "Complaint_No", "Area"), all.x=T, all.y=T)
osmgoogle$ConstituencyMatch<-ifelse(osmgoogle$OSMConstituency==osmgoogle$GoogleConstituency,1,0)
osmgoogle$OSMGoogleWardMatch<-ifelse(osmgoogle$OSMWard==osmgoogle$GoogleWard,1,0)
osmgoogle$OSMWardMatch<-ifelse(osmgoogle$OSMWard==osmgoogle$Ward,1,0)
osmgoogle$GoogleWardMatch<-ifelse(osmgoogle$GoogleWard==osmgoogle$Ward,1,0)

table(osmgoogle$ConstituencyMatch)
table(osmgoogle$OSMGoogleWardMatch)
table(osmgoogle$OSMWardMatch)
table(osmgoogle$GoogleWardMatch)
table(osmgoogle$OSMWardMatch,osmgoogle$ConstituencyMatch)

#Remove duplicates
delhicomplaints7<-merge(delhicomplaints3, osmgoogle, by=c("Ward", "Complaint_No", "Area"))
delhicomplaints8<-delhicomplaints7[!duplicated(delhicomplaints7),]
delhicomplaints8.1<-delhicomplaints8[delhicomplaints8$GoogleWardMatch==1,]

#Use Google Constituency Coding
combineddelhi<-combined[combined$City=='Delhi' & combined$Year=='2019',]
delhicomplaints9<-merge(delhicomplaints8, combineddelhi, by.x="GoogleConstituency", by.y="Constituency", all.x=T)
delhicomplaints10<-delhicomplaints9[delhicomplaints9$GoogleWardMatch==1,]
sort(table(delhicomplaints10$GoogleConstituency))


#Check representativeness of GoogleWardMatch
representation1<-delhicomplaints9[,1:34]
names(representation1)[names(representation1) == 'Ward.x'] <- 'W'
representation1<-dummy_columns(representation1, select_columns="W")


#Complaint Representativeness
category<-names(representation1[c(11:25)])
tvalues<-as.data.frame(matrix(nrow=15,ncol=5))
colnames(tvalues)<-c("name", "meanx", "meany", "tvalue", "pvalue")

for(i in 1:length(category)){
  t0<-t.test(representation1[representation1$GoogleWardMatch==1,][,i+10],
             representation1[representation1$GoogleWardMatch==0,][,i+10])
  tvalues[i,1]<-category[i]
  tvalues[i,2]<-t0$estimate[1]
  tvalues[i,3]<-t0$estimate[2]
  tvalues[i,4]<-t0$statistic
  tvalues[i,5]<-t0$p.value
}


#Ward Representativeness
ward<-names(representation1[35:46])
tvalues2<-as.data.frame(matrix(nrow=12,ncol=5))
colnames(tvalues2)<-c("name", "meanx", "meany", "tvalue", "pvalue")

for(i in 1:length(ward)){
  t0<-t.test(representation1[representation1$GoogleWardMatch==1,][,i+34],
             representation1[representation1$GoogleWardMatch==0,][,i+34])
  tvalues2[i,1]<-ward[i]
  tvalues2[i,2]<-t0$estimate[1]
  tvalues2[i,3]<-t0$estimate[2]
  tvalues2[i,4]<-t0$statistic
  tvalues2[i,5]<-t0$p.value
}


tvalues3<-as.data.frame(matrix(nrow=12,ncol=5))
colnames(tvalues3)<-c("name", "meanx", "meany", "tvalue", "pvalue")


for(i in 1:length(ward)){
  t0<-t.test(representation1[representation1$GoogleWardMatch==1,][,i+34],
             representation1[,i+34])
  tvalues3[i,1]<-ward[i]
  tvalues3[i,2]<-t0$estimate[1]
  tvalues3[i,3]<-t0$estimate[2]
  tvalues3[i,4]<-t0$statistic
  tvalues3[i,5]<-t0$p.value
}

#Table SI.3.1
print(xtable(tvalues, type='latex', row.names=F), include.rownames=FALSE)

#Table SI.3.2
print(xtable(tvalues2, type='latex', row.names=F), include.rownames=FALSE)

#Table SI.3.3
print(xtable(tvalues3, type='latex', row.names=F), include.rownames=FALSE)


#Common Area Descriptions
representation2<-as.data.frame(cbind(representation1$Area,representation1$GoogleWardMatch))
representation2$complaint<-1
representation2$V2<-as.numeric(as.character(representation2$V2))
representation3<-representation2%>%
  group_by(V1)%>%
  summarize_all(sum, na.rm=T)

#Complaint Descriptives
representation4<-as.data.frame(cbind(representation1$ComplaintCoding,representation1$GoogleWardMatch))
representation4$complaint<-1
representation4$V2<-as.numeric(as.character(representation4$V2))
representation5<-representation4%>%
  group_by(V1)%>%
  summarize_all(sum, na.rm=T)
representation5<-representation5[order(-representation5$complaint),]
colnames(representation5)<-c("ComplaintCoding", "CodedComplaints", "TotalComplaints")

#Table SI.3.4
representation6<-representation5[1:25,]
print(xtable(representation6, type='latex', row.names=F, digits=c(0)), include.rownames=FALSE)

t.test(representation5$CodedComplaints,representation5$TotalComplaints, paired=T)




#Aggregate to Constituency level
delhicomplaints11<-delhicomplaints8.1[,c(10:24,28)]
delhicomplaints11<-na.omit(delhicomplaints11)
rownames(delhicomplaints11)<-NULL
delhicomplaints11$GoogleComplaints<-1

delhicomplaints12<-delhicomplaints11%>%
  group_by(GoogleConstituency)%>%
  summarise_all(sum, na.rm=T)
colnames(delhicomplaints12)<-c("GoogleConstituency", "Complaint_Pollution","Complaint_Pests",
                                         "Complaint_Infrastructure", "Complaint_SWM",
                                         "Complaint_Roads", "Complaint_Water","Complaint_Other", "Complaint_Drainage", "Complaint_Electricity",
                                        "Complaint_Unauthorized", "Complaint_Services", "Complaint_Environment","Complaint_NA", "Complaint_Health",
                                         "Complaint_Welfare", "GoogleComplaints")
delhicomplaints12[is.na(delhicomplaints12)]<-0

combineddelhi2<-merge(combineddelhi, delhicomplaints12, by.x="Constituency", by.y="GoogleConstituency", all.x=T)


#All Questions
delhiconstituencyquestions<-delhiquestions3[,c(14,45:61)]
delhiconstituencyquestions$TotalQuestions<-1

delhiconstituencyquestions2<-delhiconstituencyquestions%>%
  group_by(Constituency)%>%
  summarise_all(sum, na.rm=T)
colnames(delhiconstituencyquestions2)<-c("Constituency", "Question_Pollution","Question_Pests",
 "Question_Infrastructure", "Question_Other", "Question_SWM", "Question_Welfare",
 "Question_Roads", "Question_Water","Question_Drainage", "Question_Education", "Question_Electricity",
 "Question_Environment", "Question_Services", "Question_Unauthorized", "Question_Health",
 "Question_NA", "Question_Revenue", "TotalQuestions")

#Remove NA row
delhiconstituencyquestions2<-delhiconstituencyquestions2[-268,]
delhiconstituencyquestions2[is.na(delhiconstituencyquestions2)]<-0

combineddelhi3<-merge(combineddelhi2, delhiconstituencyquestions2, by="Constituency", all.x=T, all.y=T)

#Only ward questions
delhiconstituencyquestions3<-delhiwardquestions3[,c(14,45:61)]
delhiconstituencyquestions3$WardQuestions<-1

delhiconstituencyquestions4<-delhiconstituencyquestions3%>%
  group_by(Constituency)%>%
  summarise_all(sum, na.rm=T)
colnames(delhiconstituencyquestions4)<-c("Constituency", "WQuestion_Pollution","WQuestion_Pests",
                                         "WQuestion_Infrastructure", "WQuestion_Other", "WQuestion_SWM", "WQuestion_Welfare",
                                         "WQuestion_Roads", "WQuestion_Water","WQuestion_Drainage", "WQuestion_Education", "WQuestion_Electricity",
                                         "WQuestion_Environment", "WQuestion_Services", "WQuestion_Unauthorized", "WQuestion_Health",
                                         "WQuestion_NA", "WQuestion_Revenue", "WardQuestions")

delhiconstituencyquestions4<-delhiconstituencyquestions4[-246,]
delhiconstituencyquestions4[is.na(delhiconstituencyquestions4)]<-0

combineddelhi4<-merge(combineddelhi3, delhiconstituencyquestions4, by="Constituency", all.x=T, all.y=T)
#Rows where there is no data (3 MCs who are mayors) are removed
combineddelhi4<-combineddelhi4[-c(34,45,259),]

#remaining NAs->0 since no questions/complaints implies 0
combineddelhi4[,c(103:154)][is.na(combineddelhi4[,c(103:154)])]<-0

combineddelhi4log<-as.data.frame(combineddelhi4[,103:154]%>%mutate_all(~log(.+0.001)))
colnames(combineddelhi4log) <- paste0(colnames(combineddelhi4log), "log")

combineddelhi4_n<-as.data.frame(combineddelhi4[,103:154]%>%mutate_all(~normalize(.)))
colnames(combineddelhi4_n) <- paste0(colnames(combineddelhi4_n), "_n")

combineddelhi4log_n<-cbind(combineddelhi4log,combineddelhi4_n)

combineddelhi5<-cbind(combineddelhi4, combineddelhi4log_n)

#save(combineddelhi5, file="combineddelhi5.RData")
#load("combineddelhi5.RData")



set.seed(123)
model0<-lm(TotalQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi5)
coef0<-coeftest(model0, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))




library(lmerTest)
model0.1<-lmerTest::lmer(TotalQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+MOV+(1|Ward),
               combineddelhi5)
summary(model0.1)

model0.2<-stan_lmer(TotalQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                      ReservedCaste+Population+C_Standing+
                      NumberCommittees+Age+BA+MOV+(1|Ward),
                    combineddelhi5)
summary(model0.2)



model1<-lm(WardQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi5)
coef1<-coeftest(model1, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model1.1<-lmerTest::lmer(WardQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+MOV+(1|Ward),
               combineddelhi5)
summary(model1.1)

model1.2<-stan_lmer(WardQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                      ReservedCaste+Population+C_Standing+
                      NumberCommittees+Age+BA+MOV+(1|Ward),
                    combineddelhi5)
summary(model1.2)

library(ggplot2)

coefs<-c(coef(summary(model0))[2,1], fixef(model0.1)[2],
         model0.2$coefficients[2], coef(summary(model1))[2,1],
         fixef(model1.1)[2], model1.2$coefficients[2])

ses<-c(coef0[2,2], sqrt(diag(vcov(model0.1, useScale = FALSE)))[2],
       model0.2$ses[2], coef1[2,2],
       sqrt(diag(vcov(model1.1, useScale = FALSE)))[2],model1.2$ses[2])
overallmodels<-as.data.frame(cbind(coefs,ses))
overallmodels$type<-factor(rep(c("Linear", "Multilevel", "Bayesian Multilevel"),2),
                           levels=c("Bayesian Multilevel", "Multilevel", "Linear"),
                           labels=c("Bayesian Multilevel", "Multilevel", "Linear"))
overallmodels$dv<-factor(c(rep("All",3),rep("Ward",3)),levels=c("Ward", "All"), labels=c("Ward", "All"))
rownames(overallmodels)<-NULL
overallmodels$coefs<-as.numeric(as.character(overallmodels$coefs))
overallmodels$ses<-as.numeric(as.character(overallmodels$ses))
overallmodels$lower<-(overallmodels$coefs-(1.96*overallmodels$ses))
overallmodels$lower[3]<-model0.2$stan_summary[2,4]
overallmodels$lower[6]<-model1.2$stan_summary[2,4]
overallmodels$upper<-(overallmodels$coefs+(1.96*overallmodels$ses))
overallmodels$upper[3]<-model0.2$stan_summary[2,10]
overallmodels$upper[6]<-model1.2$stan_summary[2,10]

overallmodels1 <- ggplot(data=overallmodels, aes(dv, coefs, color = type, 
                                                 group = type, shape = type)) + 
  geom_point(stat="identity", position = position_dodge(width = .7, 
                                                        preserve = "total"), size=3) + theme_bw()
overallmodels2 <- overallmodels1 + geom_linerange(data=overallmodels, 
                                                  aes(ymin=lower,ymax=upper), position = position_dodge(width = .7, preserve = "total"), size=1.3) + 
  labs(x="Questions",y="Complaints") + 
  coord_flip(ylim = c(-0.55, 0.55)) +  theme(axis.title=element_text(size=12)) + theme(axis.text=element_text(size=12)) + 
  theme(legend.text=element_text(size=12)) + theme(legend.title=element_text(size=12)) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme(legend.position = "right") + labs(color = "Model") + 
  labs(group = "Model") + labs(shape = "Model")

#Main text Figure 1
overallmodels2


#Table SI.5.1
library(texreg)
texreg(c(model0.1, model0.1, model0.1, model0.1), 
       override.coef=list(model0$coefficients[1:12], 
                          fixef(model0.1),
                          model1$coefficients[1:12], 
                          fixef(model1.1)),
       override.se=list(coef0[1:12,2], 
                        coef(summary(model0.1))[,2],
                        coef1[1:12,2], 
                        coef(summary(model1.1))[,2]),
       override.pvalues = list(coef0[1:12,4],
                               coef(summary(model0.1))[,4],
                               coef1[1:12,4],
                               coef(summary(model1.1))[,4]),
       stars=c(0.01, 0.05, 0.1))


#Table SI.5.2
library(xtable)
model0.2table<-as.data.frame(cbind(model0.2$coefficients[1:12], model0.2$ses[1:12]))
colnames(model0.2table)<-c("coefs", "ses")
model0.2table$lower<-model0.2$stan_summary[1:12,4]
model0.2table$upper<-model0.2$stan_summary[1:12,10]
xtable(model0.2table)

#Table SI.5.3
model1.2table<-as.data.frame(cbind(model1.2$coefficients[1:12], model1.2$ses[1:12]))
colnames(model1.2table)<-c("coefs", "ses")
model1.2table$lower<-model1.2$stan_summary[1:12,4]
model1.2table$upper<-model1.2$stan_summary[1:12,10]
xtable(model1.2table)



###Robustness check without two zero Wards
combineddelhi6<-combineddelhi5[combineddelhi5$Ward!="Keshavpuram",]
combineddelhi6<-combineddelhi6[combineddelhi6$Ward!="CivilLines",]

set.seed(123)
model0<-lm(TotalQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi6)
coef0<-coeftest(model0, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi6$Ward, type="wild"))

model0.1<-lmerTest::lmer(TotalQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+MOV+(1|Ward),
               combineddelhi6)
summary(model0.1)

model0.2<-stan_lmer(TotalQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                      ReservedCaste+Population+C_Standing+
                      NumberCommittees+Age+BA+MOV+(1|Ward),
                    combineddelhi6)
summary(model0.2)


model1<-lm(WardQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi6)
coef1<-coeftest(model1, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi6$Ward, type="wild"))


model1.1<-lmerTest::lmer(WardQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+MOV+(1|Ward),
               combineddelhi6)
summary(model1.1)

model1.2<-stan_lmer(WardQuestionslog~GoogleComplaintslog+female+SC_Pct+BJP+
                      ReservedCaste+Population+C_Standing+
                      NumberCommittees+Age+BA+MOV+(1|Ward),
                    combineddelhi6)
summary(model1.2)


#Table SI.5.4
texreg(c(model0.1, model0.1, model0.1, model0.1), 
       override.coef=list(model0$coefficients[1:12], 
                          fixef(model0.1),
                          model1$coefficients[1:12], 
                          fixef(model1.1)),
       override.se=list(coef0[1:12,2], 
                        coef(summary(model0.1))[,2],
                        coef1[1:12,2], 
                        coef(summary(model1.1))[,2]),
       override.pvalues = list(coef0[1:12,4],
                               coef(summary(model0.1))[,4],
                               coef1[1:12,4],
                               coef(summary(model1.1))[,4]),
       stars=c(0.01, 0.05, 0.1))


#Table SI.5.5
library(xtable)
model0.2table<-as.data.frame(cbind(model0.2$coefficients[1:12], model0.2$ses[1:12]))
colnames(model0.2table)<-c("coefs", "ses")
model0.2table$lower<-model0.2$stan_summary[1:12,4]
model0.2table$upper<-model0.2$stan_summary[1:12,10]
xtable(model0.2table)

#Table SI.5.6
model1.2table<-as.data.frame(cbind(model1.2$coefficients[1:12], model1.2$ses[1:12]))
colnames(model1.2table)<-c("coefs", "ses")
model1.2table$lower<-model1.2$stan_summary[1:12,4]
model1.2table$upper<-model1.2$stan_summary[1:12,10]
xtable(model1.2table)



#Individual Models for each question/complaint type
model2<-lm(Question_Pollution~Complaint_Pollution+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2<-coeftest(model2, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.1<-lm(WQuestion_Pollution~Complaint_Pollution+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.1<-coeftest(model2.1, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.2<-lm(Question_Pests~Complaint_Pests+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.2<-coeftest(model2.2, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.3<-lm(WQuestion_Pests~Complaint_Pests+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.3<-coeftest(model2.3, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.4<-lm(Question_Infrastructure~Complaint_Infrastructure+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.4<-coeftest(model2.4, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.5<-lm(WQuestion_Infrastructure~Complaint_Infrastructure+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.5<-coeftest(model2.5, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.6<-lm(Question_SWM~Complaint_SWM+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.6<-coeftest(model2.6, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.7<-lm(WQuestion_SWM~Complaint_SWM+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.7<-coeftest(model2.7, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.10<-lm(Question_Roads~Complaint_Roads+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.10<-coeftest(model2.10, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.11<-lm(WQuestion_Roads~Complaint_Roads+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.11<-coeftest(model2.11, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))



model2.12<-lm(Question_Water~Complaint_Water+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.12<-coeftest(model2.12, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.13<-lm(WQuestion_Water~Complaint_Water+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.13<-coeftest(model2.13, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))



model2.14<-lm(Question_Drainage~Complaint_Drainage+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.14<-coeftest(model2.14, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.15<-lm(WQuestion_Drainage~Complaint_Drainage+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.15<-coeftest(model2.15, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.18<-lm(Question_Electricity~Complaint_Electricity+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.18<-coeftest(model2.18, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.19<-lm(WQuestion_Electricity~Complaint_Electricity+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.19<-coeftest(model2.19, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))



model2.20<-lm(Question_Environment~Complaint_Environment+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.20<-coeftest(model2.20, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.21<-lm(WQuestion_Environment~Complaint_Environment+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.21<-coeftest(model2.21, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))



model2.22<-lm(Question_Services~Complaint_Services+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.22<-coeftest(model2.22, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model2.23<-lm(WQuestion_Services~Complaint_Services+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.23<-coeftest(model2.23, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))



model2.24<-lm(Question_Unauthorized~Complaint_Unauthorized+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.24<-coeftest(model2.24, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.25<-lm(WQuestion_Unauthorized~Complaint_Unauthorized+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.25<-coeftest(model2.25, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))



model2.26<-lm(Question_Health~Complaint_Health+female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+as.factor(Ward)
           , data=combineddelhi5)
coef2.26<-coeftest(model2.26, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))


model2.27<-lm(WQuestion_Health~Complaint_Health+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+as.factor(Ward)
             , data=combineddelhi5)
coef2.27<-coeftest(model2.27, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

#No revenue complaints


#Table SI.5.7
stargazer(model2, model2.2, model2.4, model2.6, model2.10, model2.12,
          model2.14, model2.18, model2.20, model2.22, model2.24,
          model2.26, se=list(coef2[,2], coef2.2[,2], coef2.4[,2], coef2.6[,2]
                             , coef2.10[,2], coef2.12[,2],
                             coef2.14[,2], coef2.18[,2],
                             coef2.20[,2], coef2.22[,2], coef2.24[,2], coef2.26[,2]))


#Table SI.5.8
stargazer(model2.1, model2.3, model2.5, model2.7, model2.11, model2.13,
          model2.15, model2.19, model2.21, model2.23, model2.25,
          model2.27, se=list(coef2.1[,2], coef2.3[,2], coef2.5[,2], coef2.7[,2]
                             , coef2.11[,2], coef2.13[,2],
                             coef2.15[,2], coef2.19[,2],
                             coef2.21[,2], coef2.23[,2], coef2.25[,2], coef2.27[,2]))




###Lmer models
model3<-lmer(Question_Pollution~Complaint_Pollution+female+SC_Pct+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+(1|Ward)
             , data=combineddelhi5)

model3.1<-lmer(WQuestion_Pollution~Complaint_Pollution+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.1)

model3.2<-lmer(Question_Pests~Complaint_Pests+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.2)

model3.3<-lmer(WQuestion_Pests~Complaint_Pests+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.3)


model3.4<-lmer(Question_Infrastructure~Complaint_Infrastructure+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.4)

model3.5<-lmer(WQuestion_Infrastructure~Complaint_Infrastructure+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.5)


model3.6<-lmer(Question_SWM~Complaint_SWM+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.6)

model3.7<-lmer(WQuestion_SWM~Complaint_SWM+female+SC_Pct+BJP+
                 ReservedCaste+Population+C_Standing+
                 NumberCommittees+Age+BA+
                 MOV+(1|Ward)
               , data=combineddelhi5)
summary(model3.7)

model3.10<-lmer(Question_Roads~Complaint_Roads+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.10)

model3.11<-lmer(WQuestion_Roads~Complaint_Roads+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.11)


model3.12<-lmer(Question_Water~Complaint_Water+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.12)

model3.13<-lmer(WQuestion_Water~Complaint_Water+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.13)


model3.14<-lmer(Question_Drainage~Complaint_Drainage+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.14)

model3.15<-lmer(WQuestion_Drainage~Complaint_Drainage+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.15)


model3.18<-lmer(Question_Electricity~Complaint_Electricity+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.18)

model3.19<-lmer(WQuestion_Electricity~Complaint_Electricity+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.19)


model3.20<-lmer(Question_Environment~Complaint_Environment+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.20)

model3.21<-lmer(WQuestion_Environment~Complaint_Environment+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.21)


model3.22<-lmer(Question_Services~Complaint_Services+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.22)


model3.23<-lmer(WQuestion_Services~Complaint_Services+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.23)


model3.24<-lmer(Question_Unauthorized~Complaint_Unauthorized+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.24)

model3.25<-lmer(WQuestion_Unauthorized~Complaint_Unauthorized+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.25)


model3.26<-lmer(Question_Health~Complaint_Health+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.26)

model3.27<-lmer(WQuestion_Health~Complaint_Health+female+SC_Pct+BJP+
                  ReservedCaste+Population+C_Standing+
                  NumberCommittees+Age+BA+
                  MOV+(1|Ward)
                , data=combineddelhi5)
summary(model3.27)


class(model3) <- "lmerMod"
class(model3.1) <- "lmerMod"
class(model3.2) <- "lmerMod"
class(model3.3) <- "lmerMod"
class(model3.4) <- "lmerMod"
class(model3.5) <- "lmerMod"
class(model3.6) <- "lmerMod"
class(model3.7) <- "lmerMod"
class(model3.10) <- "lmerMod"
class(model3.11) <- "lmerMod"
class(model3.12) <- "lmerMod"
class(model3.13) <- "lmerMod"
class(model3.14) <- "lmerMod"
class(model3.15) <- "lmerMod"
class(model3.18) <- "lmerMod"
class(model3.19) <- "lmerMod"
class(model3.20) <- "lmerMod"
class(model3.21) <- "lmerMod"
class(model3.22) <- "lmerMod"
class(model3.23) <- "lmerMod"
class(model3.24) <- "lmerMod"
class(model3.25) <- "lmerMod"
class(model3.26) <- "lmerMod"
class(model3.27) <- "lmerMod"

#Table SI.5.9
stargazer(model3, model3.2, model3.4, model3.6, model3.10, model3.12,
          model3.14, model3.18, model3.20, model3.22, model3.24,
          model3.26)

#Table SI.5.10
stargazer(model3.1, model3.3, model3.5, model3.7, model3.11, model3.13,
          model3.15, model3.19, model3.21, model3.23, model3.25,
          model3.27)


#Word Count
library(ngram)
wordcount(delhiquestions$Question)
620677/19537

#Look at predicting complaints
model0<-lm(GoogleComplaintslog~female+SC_Pct+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi5)
coef0<-coeftest(model0, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model0.1<-lmerTest::lmer(GoogleComplaintslog~female+SC_Pct+BJP+
                           ReservedCaste+Population+C_Standing+
                           NumberCommittees+Age+BA+MOV+(1|Ward),
                         combineddelhi5)
summary(model0.1)

model0.2<-stan_lmer(GoogleComplaintslog~female+SC_Pct+BJP+
                      ReservedCaste+Population+C_Standing+
                      NumberCommittees+Age+BA+MOV+(1|Ward),
                    combineddelhi5)
model0.2$stan_summary

#Table SI.5.11
texreg(c(model0.1, model0.1), 
       override.coef=list(model0$coefficients[1:11], 
                          fixef(model0.1)
                          ),
       override.se=list(coef0[1:11,2], 
                        coef(summary(model0.1))[,2]
                        ),
       override.pvalues = list(coef0[1:11,4],
                               coef(summary(model0.1))[,4]
                              ),
       stars=c(0.01, 0.05, 0.1))


#Table SI.5.12
library(xtable)
model0.2table<-as.data.frame(cbind(model0.2$coefficients[1:12], model0.2$ses[1:12]))
colnames(model0.2table)<-c("coefs", "ses")
model0.2table$lower<-model0.2$stan_summary[1:12,4]
model0.2table$upper<-model0.2$stan_summary[1:12,10]
xtable(model0.2table)


#Heterogeneous effects
model4<-lm(TotalQuestionslog~GoogleComplaintslog*SC_Pct+SC_Pct+female+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi5)
coef4<-coeftest(model4, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.1<-lm(TotalQuestionslog~GoogleComplaintslog*female+SC_Pct+female+BJP+
             ReservedCaste+Population+C_Standing+
             NumberCommittees+Age+BA+
             MOV+Ward
           , data=combineddelhi5)
coef4.1<-coeftest(model4.1, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.2<-lm(TotalQuestionslog~GoogleComplaintslog*BJP+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.2<-coeftest(model4.2, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.3<-lm(TotalQuestionslog~GoogleComplaintslog*ReservedCaste+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.3<-coeftest(model4.3, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.4<-lm(TotalQuestionslog~GoogleComplaintslog*Population+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.4<-coeftest(model4.4, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.5<-lm(TotalQuestionslog~GoogleComplaintslog*C_Standing+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.5<-coeftest(model4.5, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.6<-lm(TotalQuestionslog~GoogleComplaintslog*NumberCommittees+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.6<-coeftest(model4.6, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.7<-lm(TotalQuestionslog~GoogleComplaintslog*Age+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.7<-coeftest(model4.7, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.8<-lm(TotalQuestionslog~GoogleComplaintslog*BA+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.8<-coeftest(model4.8, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

model4.9<-lm(TotalQuestionslog~GoogleComplaintslog*MOV+SC_Pct+female+BJP+
               ReservedCaste+Population+C_Standing+
               NumberCommittees+Age+BA+
               MOV+Ward
             , data=combineddelhi5)
coef4.9<-coeftest(model4.9, vcov=function(x) vcovBS(x, R=5000, cluster=combineddelhi5$Ward, type="wild"))

#Table SI.5.13
stargazer(model4, model4.1, model4.2, model4.3, model4.4, model4.5, model4.6,
          model4.7, model4.8, model4.9,
          se=list(coef4[,2], coef4.1[,2], coef4.2[,2], coef4.3[,2]
                             ,coef4.4[,2], coef4.5[,2], coef4.6[,2],
                             coef4.7[,2], coef4.8[,2],
                             coef4.9[,2]))

library(interactions)
interact_female<-interact_plot(model4.1, pred="GoogleComplaintslog", modx="female",
              interval = TRUE, rug=T, legend.main="Female",
              x.label="Complaints", y.label="Questions", line.thickness = 1.5,
              int.type = "confidence", vcov=vcovBS(model4.1, R=5000, cluster=combineddelhi5$Ward, type="wild"),
              cex=1.5)+theme(text = element_text(size = 16))   

interact_population<-interact_plot(model4.4, pred="GoogleComplaintslog", modx="Population",
              interval = TRUE, rug=T, legend.main="Population",
              x.label="Complaints", y.label="Questions", line.thickness = 1.5,
              int.type = "confidence", vcov=vcovBS(model4.1, R=5000, cluster=combineddelhi5$Ward, type="wild"),
              cex=1.5)+theme(text = element_text(size = 16)) 

interact_committee<-interact_plot(model4.6, pred="GoogleComplaintslog", modx="NumberCommittees",
              interval = TRUE, rug=T, legend.main="Num. Comm.",
              x.label="Complaints", y.label="Questions", line.thickness = 1.5,
              int.type = "confidence", vcov=vcovBS(model4.1, R=5000, cluster=combineddelhi5$Ward, type="wild"),
              cex=1.5)+theme(text = element_text(size = 16)) 

interact_BA<-interact_plot(model4.8, pred="GoogleComplaintslog", modx="BA",
              interval = TRUE, rug=T, legend.main="BA",
              x.label="Complaints", y.label="Questions", line.thickness = 1.5,
              int.type = "confidence", vcov=vcovBS(model4.1, R=5000, cluster=combineddelhi5$Ward, type="wild"),
              cex=1.5)+theme(text = element_text(size = 16)) 

library(ggpubr)
#Figure SI.5.1
ggarrange(interact_female, interact_population, interact_committee, interact_BA,
          ncol = 2, nrow = 2)



#Average SC population by ward
sc<-select(combineddelhi5, c("Ward", "SC_Population", "Population"))
sc2<- sc %>%
  group_by(Ward) %>%
  summarise_all(sum)
sc2$SC_Pct<-sc2$SC_Population/sc2$Population

sink()


